ENSO analysis new

Packages

library(tidyverse)

Registered S3 methods overwritten by 'dbplyr':
  method         from
  print.tbl_lazy     
  print.tbl_sql      
── Attaching packages ───────────────────────────────────────────────────── tidyverse 1.3.2 ──✔ ggplot2 3.4.0     ✔ purrr   1.0.1
✔ tibble  3.1.8     ✔ dplyr   1.1.0
✔ tidyr   1.3.0     ✔ stringr 1.5.0
✔ readr   2.1.3     ✔ forcats 1.0.0── Conflicts ──────────────────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag()    masks stats::lag()

library(sf)

Linking to GEOS 3.11.0, GDAL 3.5.3, PROJ 9.1.0; sf_use_s2() is TRUE

library(here)

here() starts at /Users/ramkripa/Desktop/projects/emdat_proj

library(weed)
library(rnaturalearth)
library(rnaturalearthdata)

Attaching package: ‘rnaturalearthdata’

The following object is masked from ‘package:rnaturalearth’:

    countries110

Loading In the Data

emdat <- read_csv(here("full_emdat_geocoded_finally.csv"))

Warning: One or more parsing issues, call `problems()` on your data frame for details, e.g.:
  dat <- vroom(...)
  problems(dat)Rows: 74035 Columns: 51── Column specification ──────────────────────────────────────────────────────────────────────
Delimiter: ","
chr  (27): Dis No, Seq, Disaster Group, Disaster Subgroup, Disaster Type, Disaster Subtype...
dbl  (18): Year, Dis Mag Value, Start Year, Start Month, Start Day, End Year, End Month, E...
lgl   (5): Glide, Aid Contribution, Reconstruction Costs ('000 US$), Admin1 Code, uncertai...
time  (1): Local Time
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

How many disasters were located?

How many locations were located?

All types of disasters

Restrict to dry and wet hydrological disasters

dry_list <- c("Drought", "Extreme temperature", "Wildfire")
wet_list <- c("Flood", "Landslide", "Storm")

emdat_hydrological <- emdat %>%
  filter((`Disaster Type` %in% dry_list) || (`Disaster Type` %in% wet_list)) %>%
  mutate(dry = (`Disaster Type` %in% dry_list), wet = (`Disaster Type` %in% wet_list))

Warning: There were 2 warnings in `filter()`.
The first warning was:
ℹ In argument: `(`Disaster Type` %in% dry_list) || ...`.
Caused by warning in `` (`Disaster Type` %in% dry_list) || (`Disaster Type` %in% wet_list) ``:
! 'length(x) = 74035 > 1' in coercion to 'logical(1)'
ℹ Run ]8;;ide:run:dplyr::last_dplyr_warnings()dplyr::last_dplyr_warnings()]8;; to see the 1 remaining warning.

ENSO data

emdat_hydrological_enso %>%
  filter(!is.na(enso_Value)) %>%
  filter(!is.na(wet)) %>%
  mutate(enso_positive = (enso_Value > 0)) %>%
  group_by(`Dis No`) %>%
  summarize(enso_pos = any(enso_positive),enso_val = mean(enso_Value), is_wet = any(wet)) %>%
  ggplot(mapping = aes(x = enso_val, fill = is_wet)) +
    geom_histogram(stat = "density") +
    facet_grid(~is_wet) +
  ggtitle("Counts of Disasters compared to value of ONI with wet disasters in Blue and Dry in Red")

Warning: Ignoring unknown parameters: `binwidth`, `bins`, and `pad`

## calculate relative risk within each box
### create boxes

resy <- 2

emdat_hydrological_enso %>%
  left_join(event_num_months, by = c("Event")) %>%
  filter(!is.na(enso_Value)) %>%
  filter(!is.na(wet)) %>%
  filter(!is.na(lat)) %>%
  mutate(wet_factor = as_factor(wet)) %>%
  mutate(Event_factor = fct_recode(as_factor(Event), no_event = "N", negative = "-", positive = "+")) %>%
  mutate(wet_string = fct_recode(wet_factor, wet_disasters = "TRUE", dry_disasters = "FALSE")) %>%
  mutate(lat_box = (as.integer(lat) %/% resy) * resy , lng_box = (as.integer(lng) %/% resy) * resy) %>%
  mutate(dis_no_factor = as_factor(`Dis No`)) %>%
  group_by(lat_box, lng_box, Event_factor, wet_string) %>%
  summarize(num_months = first(num_months), num_disasters = n_distinct(dis_no_factor)) %>%
  mutate(disaster_rate = num_disasters / num_months) %>%
  pivot_wider(id_cols = c(lat_box, lng_box), names_from = c(Event_factor, wet_string), values_from = c(disaster_rate), values_fill = 0) %>%
  mutate(positive_wet_relative_risk = positive_wet_disasters / no_event_wet_disasters, 
         negative_wet_relative_risk = negative_wet_disasters / no_event_wet_disasters, 
         positive_dry_relative_risk =  positive_dry_disasters / no_event_dry_disasters,
         negative_dry_relative_risk = negative_dry_disasters / no_event_dry_disasters) %>%
  select(lat_box, lng_box, positive_wet_relative_risk, negative_wet_relative_risk, positive_dry_relative_risk, negative_dry_relative_risk) %>%
  pivot_longer(cols = c(positive_wet_relative_risk, negative_wet_relative_risk, positive_dry_relative_risk, negative_dry_relative_risk), names_to = 'risk', values_to = 'value') %>%
  mutate(risk_factor = as_factor(risk), values_cleaned = if_else(is.infinite(value), 18, value)) %>%
  mutate(relative_risk_cleaned = if_else(is.nan(values_cleaned), 0, values_cleaned)) %>%
  ggplot() +
  geom_sf(data = world_map) +
  geom_tile(mapping = aes(x = lng_box, y = lat_box, fill = relative_risk_cleaned, width = resy)) +
  scale_fill_gradient(low = 'white', high = 'red') +
  facet_grid(rows = vars(risk_factor))+
  ggtitle("Relative risk plot using calculations of the form [{num_disasters_positive_wet / num_months_positive} / {num_disasters_no_event_wet / num_months_no_event}]")

`summarise()` has grouped output by 'lat_box', 'lng_box', 'Event_factor'. You can override using the `.groups` argument.

## calculate percent of total disasters within each box
### create boxes 2x2

resy <- 2

emdat_hydrological_enso %>%
  left_join(event_num_months, by = c("Event")) %>%
  filter(!is.na(enso_Value)) %>%
  filter(!is.na(wet)) %>%
  filter(!is.na(lat)) %>%
  mutate(wet_factor = as_factor(wet)) %>%
  mutate(Event_factor = fct_recode(as_factor(Event), no_event = "N", negative = "-", positive = "+")) %>%
  mutate(wet_string = fct_recode(wet_factor, wet_disasters = "TRUE", dry_disasters = "FALSE")) %>%
  mutate(lat_box = (as.integer(lat) %/% resy) * resy , lng_box = (as.integer(lng) %/% resy) * resy) %>%
  mutate(dis_no_factor = as_factor(`Dis No`)) %>%
  group_by(lat_box, lng_box, Event_factor, wet_string) %>%
  summarize(num_months = first(num_months), num_disasters = n_distinct(dis_no_factor)) %>%
  mutate(disaster_rate = num_disasters / num_months) %>%
  pivot_wider(id_cols = c(lat_box, lng_box), names_from = c(Event_factor, wet_string), values_from = c(disaster_rate), values_fill = 0) %>%
  mutate(positive_wet_percent = positive_wet_disasters / (no_event_wet_disasters + positive_wet_disasters + negative_wet_disasters) , 
         negative_wet_percent = negative_wet_disasters / (no_event_wet_disasters + positive_wet_disasters + negative_wet_disasters), 
         positive_dry_percent =  positive_dry_disasters / (no_event_dry_disasters + positive_dry_disasters + negative_dry_disasters),
         negative_dry_percent =  negative_dry_disasters / (no_event_dry_disasters + positive_dry_disasters + negative_dry_disasters)
         ) %>%
  select(lat_box, lng_box, positive_wet_percent, negative_wet_percent, positive_dry_percent, negative_dry_percent) %>%
  pivot_longer(cols = c(positive_wet_percent, negative_wet_percent, positive_dry_percent, negative_dry_percent), names_to = 'typey', values_to = 'value') %>%
  mutate(type_factor = as_factor(typey)) %>%
  ggplot() +
  geom_sf(data = world_map) +
  geom_tile(mapping = aes(x = lng_box, y = lat_box, fill = value, width = resy)) +
  scale_fill_gradient(low = 'white', high = 'red') +
  facet_grid(rows = vars(type_factor))+
  ggtitle("Disaster Percent plot using calculations of the form [positive_wet / all_wet]")

`summarise()` has grouped output by 'lat_box', 'lng_box', 'Event_factor'. You can override using the `.groups` argument.